!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000, 2003 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!###########################################################################

SUBROUTINE timestep()

  USE mem_basic, ONLY: &
       basic_g  ! INTENT(INOUT)

  USE node_mod, ONLY:  &
       mzp, mxp, myp,  & ! INTENT(IN)
       ia, iz, ja, jz, & ! INTENT(IN)
       i0, j0,         & ! INTENT(IN)
       izu, jzv,       & ! INTENT(IN)
       mynum,          & ! INTENT(IN)
       ibcon,          & ! INTENT(IN)
       ipara,          & ! INTENT(IN)
       mmxp,           & ! (IN)
       mmyp,           & ! (IN)
       maxgrds           ! (IN)

  !use mem_radiate, only: ! Not used

  USE mem_cuparm, ONLY: &
       NNQPARM, & ! INTENT(IN)
       IF_CUINV   ! INTENT(IN)

  USE mem_varinit, ONLY: &
       NUD_TYPE ! INTENT(IN)

  USE mem_turb, ONLY: &
       IF_URBAN_CANOPY, & ! INTENT(IN)
       ihorgrad           ! INTENT(IN)

  USE mem_oda,   ONLY: &
       if_oda ! INTENT(IN)

  USE micphys,   ONLY: &
       level ! INTENT(IN)

  USE mem_grid, ONLY: &
       ngrid,      & ! INTENT(IN)
       time,       & ! INTENT(IN)
       dtlong,     & ! INTENT(IN)
       dtlongn,    & ! INTENT(IN)
       iyear1,     & ! INTENT(IN)
       imonth1,    & ! INTENT(IN)
       idate1,     & ! INTENT(IN)
       grid_g,     & ! INTENT(INOUT)
       nxtnest,    & ! INTENT(IN)
       if_adap,    & ! INTENT(IN)
       dtlt,       & ! INTENT(IN)
       istp,       & ! INTENT(IN)
       jdim,       & ! INTENT(IN)
       nzp,        & ! INTENT(IN)
       f_thermo_e, & ! INTENT(IN)
       f_thermo_w, & ! INTENT(IN)
       f_thermo_s, & ! INTENT(IN)
       f_thermo_n, & ! INTENT(IN)
       itime1,     &
       zt,         &
       zm,         &
       dzt,        &
       ngrids

  USE shcu_vars_const, ONLY: & ! For Shallow Cumulus Paramet.
       NNSHCU ! INTENT(IN)

  USE mem_scalar, ONLY: & ! For SiB
       scalar_g ! INTENT(IN)

  USE mem_leaf, ONLY: & ! For SiB
       ISFCL ! INTENT(IN)

  ! CAT
  USE catt_start, ONLY: &
       CATT ! INTENT(IN)

!--(DMK-CCATT)---------------------------------------------------------
!--(DMK-original)------------------------------------------------------
!  ! For CATT
!  !!use burns, only : queimadas
!  use emission_source_map, only: &
!       burns ! Subroutine
!--(DMK-CCATT-END)-----------------------------------------------------

 ! TEB_SPM
  USE teb_spm_start, ONLY: &
       TEB_SPM ! INTENT(IN)

  ! For TEB_SPM
  USE mem_emiss, ONLY: &
       ichemi,         & ! INTENT(IN)
       isource           ! INTENT(IN)

  ! ALF
  ! Necessary in new advection scheme
  USE advect_kit, ONLY : &
       calc_advec          ! Subroutine

  ! For specific optimization depending the type of machine
  USE machine_arq, ONLY: &
       machine ! INTENT(IN)

!--(DMK-CCATT)---------------------------------------------------------
  !-srf- stilt
  USE mem_stilt, ONLY: &
       stilt_g,        &
       iexev,          &
       imassflx

  !-srf for dry dep
  USE  module_dry_dep, ONLY: &
       drydep_driver           ! Subroutine
  
  !-srf - chem
  USE module_chemistry_driver, ONLY: &
       chemistry_driver        ! Subroutine

  USE chem_sources, ONLY :  &
       sources_driver          ! Subroutine

  USE mem_chem1, ONLY: &
       nvert_src=>chem1_src_z_dim_g, &
       chem1_g,                      &
       nsrc,                         & ! (IN)
       bburn,                        & ! (IN)
       chem1_src_g,                  & ! %sc_src(INOUT)
       antro,                        & ! (IN)
       bioge,                        & ! (IN)
       CHEMISTRY,                    &
       N_DYN_CHEM,                   &
       split_method,                 &
       ntimes_src,                   &
       diur_cycle

  USE mem_radiate, ONLY: &
       radiate_g

  USE grid_dims, ONLY: &
       nzpmax            ! (IN)

  USE rconstants , ONLY: &
       pi180            ! (IN)

  USE plumegen_coms,ONLY:   &
       nkp,                 & ! (IN)
       ntime,               & ! (IN)
       zt_plumegen=>zt,     & ! (INOUT)
       zm_plumegen=>zm,     & ! (INOUT)
       dzt_plumegen=>dzt,   & ! (INOUT)
       dzm,                 & ! (INOUT)
       dz,                  & ! (INOUT)
       ncall,               & ! (INOUT)
       nrec,                & ! (INOUT)
       n_setgrid              ! (INOUT)

  USE chem1_list, ONLY :          &
       chem_nspecies =>nspecies,  & ! (IN)
       spc_chem_alloc=>spc_alloc, & ! (IN)
       spc_chem_name =>spc_name,  & ! (IN)
       src,                       & ! (IN)
       on,                        & ! (IN)
       off,                       & ! (IN)
       transport

  USE aer1_list, ONLY :          &
       aer_nspecies=>nspecies,   & ! (IN)
       spc_aer_alloc=>spc_alloc, & ! (IN)
       nmodes,                   & ! (IN)
       aer_bburn => bburn,       & ! (IN)
       aer_sdust => sdust,       &
       aer_urban => urban,       &
       aer_bioge => bioge,       &
       aer_marin => marin 

  USE mem_aer1, ONLY:                   &
       AEROSOL,                         & ! (IN)
       aer1_g,                          & ! %sc_src(INOUT)
       aer_nvert_src=>aer1_src_z_dim_g

  USE mem_plume_chem1, ONLY: &
       nveg_agreg,           & ! (IN)
       tropical_forest,      & ! (IN)
       boreal_forest,        & ! (IN)
       savannah,             & ! (IN)
       grassland,            & ! (IN)
       plume_mean_g            ! %flam_frac(IN), %fire_size(IN)

!--(DMK-CCATT-END)-----------------------------------------------------

  IMPLICIT NONE

  REAL :: t1,w1
  REAL, EXTERNAL :: cputime

!  integer, save :: ncall=0

  INTEGER, PARAMETER :: acct = 0 ! To Activate acctimes

!--(DMK-CCATT)---------------------------------------------------------
  !- para testar acoplamento shallow/deep 
  INTEGER, PARAMETER :: acoshdp = 0 
!--(DMK-CCATT-END)-----------------------------------------------------

  IF (acct/=0) CALL acctimes('init',0,' ',t1,w1)

  !        +-------------------------------------------------------------+
  !        |   Timestep driver for the hybrid non-hydrostatic time-split |
  !        |      model.                                                 |
  !        +-------------------------------------------------------------+

  !  Zero out all tendency arrays.   
  !--------------------------------
  t1=cputime(w1)
  CALL TEND0()          
  IF (acct/=0) CALL acctimes('accu',1,'TEND0',t1,w1)


!--(DMK-CCATT)---------------------------------------------------------
  ! CATT & Chemistry == CCATT
  !----------------------------------------
  IF (CATT==1) THEN

     t1=cputime(w1)
     !- emission 

!     call sources_driver(mzp,mxp,myp,ia,iz,ja,jz)

     CALL sources_driver(mzp,mxp,myp,ia,iz,ja,jz, &
                         ngrid,mmxp,mmyp,maxgrds,radiate_g(ngrid)%cosz,basic_g(ngrid)%theta, &
                         basic_g(ngrid)%pp,basic_g(ngrid)%pi0,basic_g(ngrid)%rv,             &
                         basic_g(ngrid)%dn0,time,iyear1,imonth1,idate1,itime1,dtlt,          &
                         grid_g(ngrid)%rtgt,grid_g(ngrid)%lpw,grid_g(ngrid)%glat,            &
                         grid_g(ngrid)%glon,zt,zm,dzt,ngrids,nzpmax,pi180,nkp,ntime,         &
                         zt_plumegen,zm_plumegen,dzt_plumegen,dzm,dz,ncall,nrec,n_setgrid,   &
                         chem_nspecies,spc_chem_alloc,spc_chem_name,src,on,off,transport,    &
                         aer_nspecies,spc_aer_alloc,nmodes,aer_bburn,aer_sdust,aer_urban,    &
                         aer_bioge,aer_marin,nvert_src,chem1_g,nsrc,bburn,chem1_src_g,antro, &
                         bioge,chemistry,ntimes_src,diur_cycle,aerosol,aer1_g,aer_nvert_src, &
                         nveg_agreg,tropical_forest,boreal_forest,savannah,grassland,        &
                         plume_mean_g,stilt_g(ngrid)%dnp,iexev)

     IF (acct/=0) CALL acctimes('accu',1,'SOURCES',t1,w1)     
     
     !- large and subgrid scale forcing for shallow and deep cumulus
     IF( NNQPARM(ngrid).EQ.2 .AND. NNSHCU(ngrid)==2 ) CALL prepare_lsf(1)

  ENDIF
!--(DMK-CCATT-END)-----------------------------------------------------

  !  Thermodynamic diagnosis   
  !--------------------------------
  t1=cputime(w1)
  IF (level  /=  3) THEN
     CALL THERMO(mzp,mxp,myp,ia,iz,ja,jz,'SUPSAT') 
  ENDIF
  IF (acct/=0) CALL acctimes('accu',2,'THERMO',t1,w1)

!--(DMK-CCATT)---------------------------------------------------------
  !ML/SRF---------- David's mass conservation fix ------------------------
  !- Right before calling "RADIATE", add a function call:
  IF(iexev == 2)CALL exevolve(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,dtlt,'ADV')
  !-------------------------------- 
!--(DMK-CCATT-END)-----------------------------------------------------


  !  Radiation parameterization
  !--------------------------------
  t1=cputime(w1)
  CALL RADIATE(mzp,mxp,myp,ia,iz,ja,jz,mynum) 
  IF (acct/=0) CALL acctimes('accu',3,'RADIATE',t1,w1)

  !  Surface layer, soil and veggie model
  !----------------------------------------
  t1=cputime(w1)
  IF (isfcl <= 2) THEN
     CALL SFCLYR(mzp,mxp,myp,ia,iz,ja,jz,ibcon)
  ELSEIF (isfcl == 3) THEN
     CALL sfclyr_sib(mzp,mxp,myp,ia,iz,ja,jz,ibcon)
  ENDIF
  IF (acct/=0) CALL acctimes('accu',4,'SFCLYR',t1,w1)

  ! SiB for BRAMS
  ! CO2  bio source
  IF (ISFCL == 3) THEN
     t1=cputime(w1)
     CALL co2_biosource(mzp,mxp,myp,ia,iz,ja,jz,ngrid,  &
          scalar_g(1,ngrid)%sclt(1),basic_g(ngrid)%dn0,grid_g(ngrid)%rtgt)
     IF (acct/=0) CALL acctimes('accu',4,'co2_biosource',t1,w1)
  ENDIF
  !----------------------------------------
  !- CATT & Chemistry == CCATT
  !----------------------------------------
   IF (CATT==1) THEN     
     !srf- dry deposition of gases/particles
     CALL drydep_driver(mzp,mxp,myp,ia,iz,ja,jz)     
     
     !srf- large and subgrid scale forcing for shallow and deep cumulus
     IF( NNQPARM(ngrid).EQ.2 .AND. NNSHCU(ngrid)==2 ) CALL prepare_lsf(2)

  ENDIF

  !  Send boundaries to adjoining nodes
  !-------------------------------------------
  t1=cputime(w1)
  IF (ipara  ==  1) THEN
     CALL node_sendlbc() !!_f2c()
     IF (ngrid  ==  1) CALL node_sendcyclic(1) !!_f2c(1)
  ENDIF
  IF (acct/=0) CALL acctimes('accu',26,'SendLBC',t1,w1)

  !  Coriolis terms
  !  ----------------------------------------
  t1=cputime(w1)
  CALL CORLOS(mzp,mxp,myp,i0,j0,ia,iz,ja,jz,izu,jzv) 
  IF (acct/=0) CALL acctimes('accu',5,'CORLOS',t1,w1)

  !  Velocity advection
  !----------------------------------------
  t1=cputime(w1)
  ! Use Optmized advection only in SX-6, for the moment
  IF (machine==0) THEN
     ! If Generic IA32 use old Advction Scheme
     CALL ADVECTc('V',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum)
  ELSEIF (machine==1) THEN
     ! Using optmized advection scheme only in SX-6
     CALL calc_advec('V',ngrid,mzp,mxp,myp)
  ENDIF
  IF (acct/=0) CALL acctimes('accu',6,'ADVECTv',t1,w1)

  !  Cumulus parameterization
  !----------------------------------------
  t1=cputime(w1)
  IF(NNQPARM(ngrid) == 1 .OR. IF_CUINV == 1) CALL cuparm()      
  IF (acct/=0) CALL acctimes('accu',7,'CUPARM',t1,w1)

  !  Urban canopy parameterization
  !----------------------------------------
  t1=cputime(w1)
  IF( IF_URBAN_CANOPY == 1) CALL urban_canopy()      
  IF (acct/=0) CALL acctimes('accu',7,'CUPARM',t1,w1)

  !  Analysis nudging and boundary condition
  !------------------------------------------
  t1=cputime(w1)
  IF(NUD_TYPE > 0) CALL DATASSIM()  
  IF (acct/=0) CALL acctimes('accu',8,'DATASSIM',t1,w1)

  !  Observation data assimilation 
  !----------------------------------------
  t1=cputime(w1)
  IF(IF_ODA == 1) CALL oda_nudge()  
  IF (acct/=0) CALL acctimes('accu',8,'DATASSIM',t1,w1)

  !  Nested grid boundaries
  !----------------------------------------
  t1=cputime(w1)
  IF(nxtnest(ngrid).GE.1) CALL nstbdriv()  
  IF (acct/=0) CALL acctimes('accu',9,'NSTBDRIV',t1,w1)

  !  Rayleigh friction for theta
  !----------------------------------------
  t1=cputime(w1)
  CALL RAYFT()           
  IF (acct/=0) CALL acctimes('accu',11,'RAYFT',t1,w1)

  !  Get the overlap region between parallel nodes
  !---------------------------------------------------
  t1=cputime(w1)
  IF(ipara == 1) THEN      
     CALL node_getlbc() !!_f2c()  
     IF (ngrid  ==  1) CALL node_getcyclic(1) !!_f2c(1)
  ENDIF
  IF (acct/=0) CALL acctimes('accu',13,'GETlbc',t1,w1)

  !  Sub-grid diffusion terms
  !----------------------------------------
  t1=cputime(w1)
  IF ((if_adap==0).AND.(ihorgrad==2)) THEN
     CALL diffuse_brams31() !call optimized subroutine
  ELSE
     CALL diffuse()
  ENDIF
  IF (acct/=0) CALL acctimes('accu',12,'DIFFUSE',t1,w1)

  !- CATT & Chemistry == CCATT
  !----------------------------------------
  IF (CATT==1) THEN
     !- large and subgrid scale forcing for shallow and deep cumulus
     !- neste ponto LSF = radiation + pbl(vertical) diffusion,
     !s- este eh o forcing para o shallow
     IF( NNQPARM(ngrid).EQ.2 .AND. NNSHCU(ngrid)==2 ) CALL prepare_lsf(3)
     
  ENDIF

  !  Velocity advection
  !----------------------------------------
  t1=cputime(w1)
  ! Use Optmized advection only in SX-6, for the moment
  IF (machine==0) THEN
     ! If Generic IA32 use old Advction Scheme
     CALL ADVECTc('T',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum)
  ELSEIF (machine==1) THEN
     ! Using optmized advection scheme only in SX-6
     IF (ngrid<=2) THEN
        CALL calc_advec('T',ngrid,mzp,mxp,myp)
     ELSE
        CALL ADVECTc('T',mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum)
     ENDIF
  ENDIF
  IF (acct/=0) CALL acctimes('accu',15,'ADVECTs',t1,w1)

!--(DMK-CCATT)---------------------------------------------------------
  ! STILT-BRAMS coupling (ML)
  IF (imassflx == 1) CALL prep_advflx_to_stilt(mzp,mxp,myp,ia,iz,ja,jz,ngrid)
!--(DMK-CCATT-END)-----------------------------------------------------

  IF (CATT==1) THEN

!--(DMK-CCATT)---------------------------------------------------------
     !- large and subgrid scale forcing for shallow and deep cumulus
     IF( NNQPARM(ngrid).EQ.2 .AND. NNSHCU(ngrid)==2 ) CALL prepare_lsf(4)

     IF(acoshdp == 0 ) THEN

     !- cumulus parameterization by G. Grell 
     !                    Deep Convection scheme
     !- call deep first, if there is deep convection , turn off shallow.
          IF(NNQPARM(ngrid)==2) CALL CUPARM_GRELL_CATT(1) 
     !
     !                    Shallow Convection scheme
         IF(NNSHCU(ngrid)==2 .AND. NNQPARM(ngrid).EQ.2) CALL CUPARM_GRELL_CATT(2) 

     ELSEIF(acoshdp == 1 ) THEN

     !                    Shallow Convection scheme
        IF(NNSHCU(ngrid)==2 .AND. NNQPARM(ngrid).EQ.2) CALL CUPARM_GRELL_CATT(2) 

        CALL prepare_lsf(5)
        
	!- call deep first, if there is deep convection , turn off shallow.
        IF(NNQPARM(ngrid)==2) CALL CUPARM_GRELL_CATT(1) 

     ENDIF

     !- task 2:  NO production by "eclair" 
     CALL chemistry_driver(mzp,mxp,myp,ia,iz,ja,jz,2,50)
  ENDIF

!--(DMK-CCATT)---------------------------------------------------------
  !- CATT & Chemistry == CCATT
  !----------------------------------------
  IF (CATT==1 .AND. split_method== 'PARALLEL' .AND. N_DYN_CHEM == 1) THEN
     t1=cputime(w1)
     ! task 3 : production/loss by chemical processes and inclusion of the 
     ! chemistry tendency at the total tendency
     CALL chemistry_driver(mzp,mxp,myp,ia,iz,ja,jz,3,50)
     IF (acct/=0) CALL acctimes('accu',1,'CHEMISTRY',t1,w1)
  ENDIF
  IF (CATT==1 ) THEN
  !- MP 21/02/08
     t1=cputime(w1)
     ! task 4 : mass transfer between gas and liquid
     call chemistry_driver(mzp,mxp,myp,ia,iz,ja,jz,4,50)
     IF (acct/=0) CALL acctimes('accu',1,'CHEMISTRY',t1,w1)
  ENDIF
!     !---------------------------------------------------
!  ! Shallow  cumulus parameterization by Souza
   t1=cputime(w1)
   IF(NNSHCU(ngrid) == 1) CALL SHCUPA()
   CALL acctimes(19,'SHCUPARM',t1,w1)
   IF (acct/=0) CALL acctimes('accu',26,'SHCUPARM',t1,w1)
!  !---------------------------------------------------

  IF (TEB_SPM==1) THEN
     ! Update urban emissions
     !----------------------------------------
     IF(isource==1)THEN
        t1=cputime(w1)
        CALL sources_teb(mzp,mxp,myp,ia,iz,ja,jz,ngrid,dtlt)
        CALL acctimes('accu',28,'EMISS',t1,w1)
     ENDIF
     !  Update chemistry
     !----------------------------------------
     IF(ichemi==1)THEN
        t1=cputime(w1)
        CALL ozone(mzp,mxp,myp,ia,iz,ja,jz,ngrid,dtlt)
        CALL acctimes('accu',29,'OZONE',t1,w1)
     ENDIF
  ENDIF

  !  Update scalars
  !----------------------------------------
  t1=cputime(w1)
  CALL PREDTR()          
  IF (acct/=0) CALL acctimes('accu',16,'PREDTR',t1,w1)

  !  Moisture variables positive definite
  !----------------------------------------
  t1=cputime(w1)
  CALL negadj1(mzp,mxp,myp) 
  IF (acct/=0) CALL acctimes('accu',17,'NEGADJ1',t1,w1)

  !  Microphysics
  !----------------------------------------
  t1=cputime(w1)
  IF (level==3) THEN
     IF (machine==1 .AND. TEB_SPM==0) THEN
        ! Optimized version only for SX-6
        CALL micro_opt()
     ELSE
        ! Original Version used in a Generic IA32 machine
        CALL micro()
     ENDIF
  ENDIF
  IF (acct/=0) CALL acctimes('accu',19,'MICRO',t1,w1)


 !MP-0808
 !- chemistry - microphysics tranfers - sedimentation and tranfer from clouds to rain
  !----------------------------------------
  IF (CATT==1) THEN
     t1=cputime(w1)
     ! task 5 : sedimentation and mass transfer between clouds and rain 
     call chemistry_driver(mzp,mxp,myp,ia,iz,ja,jz,5,50)
     IF (acct/=0) CALL acctimes('accu',1,'CHEMISTRY',t1,w1)
  ENDIF
 ! end change MP_0808     


  !  Thermodynamic diagnosis
  !----------------------------------------
  t1=cputime(w1)
  IF (level /= 3) THEN
     CALL THERMO(mzp,mxp,myp,1,mxp,1,myp,'MICRO') 
  ENDIF
  IF (acct/=0) CALL acctimes('accu',20,'THERMO',t1,w1)

!--(DMK-CCATT)---------------------------------------------------------
  !ML/SRF -> David's mass conservation fix ------
  ! Right before calling "TRSETS", add a function call:
  IF(iexev == 2)CALL exevolve(mzp,mxp,myp,ia,iz,ja,jz,izu,jzv,mynum,dtlt,'THV')
!--(DMK-CCATT-END)-----------------------------------------------------

  !  Apply scalar b.c.'s
  !----------------------------------------
  t1=cputime(w1)
  CALL TRSETS()          
  IF (acct/=0) CALL acctimes('accu',21,'TRSETS',t1,w1)

  !  Lateral velocity boundaries - radiative
  !-------------------------------------------
  t1=cputime(w1)
  CALL LATBND()
  IF (acct/=0) CALL acctimes('accu',10,'LATBND',t1,w1)

  !  First stage Asselin filter
  !----------------------------------------
  t1=cputime(w1)
  CALL HADVANCE(1)     
  IF (acct/=0) CALL acctimes('accu',22,'HADVANCE',t1,w1)

  !  Buoyancy term for w equation
  !----------------------------------------
  t1=cputime(w1)
  CALL BUOYANCY()
  IF (acct/=0) CALL acctimes('accu',23,'BUOYANCY',t1,w1)

  !  Acoustic small timesteps
  !----------------------------------------
  t1=cputime(w1)
  !call ACOUSTIC()
  CALL ACOUSTIC_new()
  IF (acct/=0) CALL acctimes('accu',24,'ACOUSTIC',t1,w1)

  !  Last stage of Asselin filter
  !----------------------------------------
  t1=cputime(w1)
  CALL HADVANCE(2)

  !  Velocity/pressure boundary conditions
  !----------------------------------------
  CALL VPSETS()          
  IF (acct/=0) CALL acctimes('accu',25,'HADVANCE',t1,w1)

  IF (acct/=0) THEN
     IF(MOD(istp,4) == 0) THEN
        CALL acctimes('prin',-1,' ',t1,w1)
     ENDIF
  ENDIF

!--(DMK-CCATT)---------------------------------------------------------
  !srf-chem:  get the true air density for chemistry 
  !           to assure mass conservation 
  IF(iexev == 2) CALL get_true_air_density(mzp,mxp,myp,ia,iz,ja,jz)
  !----------------------------------------
!--(DMK-CCATT-END)-----------------------------------------------------

  ! Call THERMO on the boundaries
  CALL thermo_boundary_driver((time+dtlongn(ngrid)), dtlong, &
       f_thermo_e, f_thermo_w, f_thermo_s, f_thermo_n, &
       nzp, mxp, myp, jdim)



 !- CATT & Chemistry == CCATT
  !----------------------------------------
  IF (CATT==1) THEN 
    IF( (split_method== 'PARALLEL' .AND. N_DYN_CHEM > 1) .OR.  &
        (split_method== 'SEQUENTIAL'                   ) .OR.  &
	(split_method== 'SYMMETRIC'                    )       )THEN

       t1=cputime(w1)

       ! task 3 : production/loss by chemical processes and final updated
       !  	of each specie
       CALL chemistry_driver(mzp,mxp,myp,ia,iz,ja,jz,3,50)
       IF (acct/=0) CALL acctimes('accu',1,'CHEMISTRY',t1,w1)
    ENDIF
 ENDIF



!!$  call mass_flux(nzp,nxp,nyp,mzp,mxp,myp,a(iup),a(ivp),a(iwp) &
!!$       ,a(idn0),a(irtgu),a(irtgv),a(idyu),a(idxv),a(ipp),a(ipi0))

  RETURN
END SUBROUTINE timestep

!*************************************************************************

SUBROUTINE acctimes(action,num,string,t1,w1)

  USE mem_all
  USE node_mod

  IMPLICIT NONE
  INTEGER :: num
  REAL :: t1,w1
  CHARACTER(len=*) :: string,action
  INTEGER, PARAMETER :: num_times=100
  REAL, SAVE :: rtimes(num_times),wtimes(num_times)
  CHARACTER(len=8),SAVE :: crtimes(num_times)

  REAL, EXTERNAL :: cputime,walltime,valugp
  REAL :: sumtime,pcpu,fsecs,cpuinc,ww
  INTEGER :: i,j,npts,ip,jp,kp

  IF(action(1:4)=='init') THEN
     crtimes(1:num_times)=' '
     rtimes(1:num_times)=0.
     wtimes(1:num_times)=0.
     basic_g(ngrid)%cputime(1:mxp,1:myp)=0.
  ELSEIF(action(1:4)=='prin') THEN
     IF(mynum/=2.AND.mynum/=0) RETURN
     sumtime=0.
     DO i=1,num_times
        sumtime=sumtime+rtimes(i)
     ENDDO

     WRITE(6,*) '======= total CPU =====',sumtime,'=========='
     DO i=1,26
        WRITE(6,'(a10,i4,''-'',a12,f10.3,f7.2,2f9.3)') &
             'Timings-',mynum,crtimes(i)  &
             ,rtimes(i),rtimes(i)/sumtime*100.,wtimes(i)  &
             ,wtimes(i)-rtimes(i)
     ENDDO
     sumtime=0.
     DO j=1,myp
        DO i=1,mxp
           sumtime=sumtime+basic_g(ngrid)%cputime(i,j)
        ENDDO
     ENDDO
     WRITE(6,'(a,2i5,f10.5)') 'Total CPU secs -ngrid,mynum,secs:'  &
          ,ngrid,mynum,sumtime

  ELSEIF(action(1:4)=='accu'.OR.action(1:4)=='null') THEN
     ! Accumulate full times into tables
     pcpu=(cputime(ww)-t1)
     rtimes(num)=rtimes(num) + pcpu
     crtimes(num)=string

     fsecs=72559200.
!     wtimes(num)=wtimes(num)+(walltime(fsecs)-w1)
     wtimes(num)=wtimes(num)+(ww-w1)

     IF(action(1:4)=='accu') THEN
        ! Divide cpu time equally amoung columns and accumulate in
        ! basic_g(ngrid)%cputime(1,1)
        npts=(iz-ia+1)*(jz-ja+1)
        cpuinc=pcpu/float(npts)
        DO j=ja,jz
           DO i=ia,iz
              basic_g(ngrid)%cputime(i,j)=basic_g(ngrid)%cputime(i,j)+cpuinc
           ENDDO
        ENDDO
     ENDIF
  ENDIF

  !    only here for debugging purposes
  ip=2
  jp=2
  kp=2
  !   if( (mynum == 1.or.mynum == 0).and.ngrid == 2) then
  ! if( (mynum == 1.or.mynum == 0)) then
  !    print '(a10,2i3,f9.1,12e14.7)',string,2,mynum,time  &
       !       ,radiate_g(ngrid)%rshort(ip,jp) &
  !       ,radiate_g(ngrid)%rlongup(ip,jp) 
  !       ,basic_g(ngrid)%up(kp,ip,jp)  &
       !       ,valugp(mzp,mxp,myp,kp,ip,jp,tend%ut(1))  &
  !       ,basic_g(ngrid)%vp(kp,ip,jp)  &
       !       ,valugp(mzp,mxp,myp,kp,ip,jp,tend%vt(1))  &
  !       ,basic_g(ngrid)%wp(kp,ip,jp)  &
       !       ,valugp(mzp,mxp,myp,kp,ip,jp,tend%wt(1)) &
  !       ,basic_g(ngrid)%pp(kp,ip,jp)  &
       !       ,valugp(mzp,mxp,myp,kp,ip,jp,tend%pt(1))  &
  !       ,basic_g(ngrid)%pi0(kp,ip,jp)  
  !       ,basic_g(ngrid)%theta(kp,ip,jp) &
       !       ,basic_g(ngrid)%thp(kp,ip,jp) 
  !       ,float(grid_g(ngrid)%lpw(ip,jp))  &
       !              ,basic_g(2)%rtp(kp,ip,jp)  &
  !              ,basic_g(2)%rv(kp,ip,jp)  &
  !              ,valugp(mzp,mxp,myp,kp,ip,jp,tend%rtt(1)) 

  !  do jp=22,1,-1
  !     print'(i3,20g12.4)',jp,(basic_g(ngrid)%pp(kp,ip,jp)),ip=1,8)
  !  enddo

  !  endif
  RETURN
END SUBROUTINE acctimes

!*************************************************************************

SUBROUTINE mass_flux(n1,n2,n3,m1,m2,m3,up,vp,wp  &
     ,dn0,rtgu,rtgv,dyu,dxv,pp,pi0)

  USE mem_grid
  USE rconstants

  IMPLICIT NONE
  INTEGER :: n1,n2,n3,m1,m2,m3
  REAL :: up(m1,m2,m3),vp(m1,m2,m3),wp(m1,m2,m3)  &
       ,dn0(n1,n2,n3),rtgu(n2,n3),dyu(n2,n3),dxv(n2,n3)  &
       ,rtgv(n2,n3),pp(m1,m2,m3),pi0(n1,n2,n3)

  REAL, SAVE :: aintmass=0.

  INTEGER :: i,j,k
  REAL :: wmass,emass,smass,nmass,prtot,tmass,ppp,area

  !cc      if (mod(time,300.).gt..1) return

  !  west/east bound
  wmass=0.
  emass=0.
  DO j=2,nyp-1
     DO k=2,nzp-1
        i=1
        wmass=wmass +  &
             up(k,i,j)*rtgu(i,j)/(dyu(i,j)*dzt(k))  &
             *(dn0(k,i,j)+dn0(k,i+1,j))*.5
        i=nxp-1
        emass=emass -  &
             up(k,i,j)*rtgu(i,j)/(dyu(i,j)*dzt(k))  &
             *(dn0(k,i,j)+dn0(k,i+1,j))*.5
     ENDDO
  ENDDO

  !  north/south bound
  smass=0.
  nmass=0.
  DO i=2,nxp-1
     DO k=2,nzp-1
        j=1
        smass=smass +  &
             vp(k,i,j)*rtgv(i,j)/(dxv(i,j)*dzt(k))  &
             *(dn0(k,i,j)+dn0(k,i,j+1))*.5
        j=nyp-1
        nmass=nmass -  &
             vp(k,i,j)*rtgv(i,j)/(dxv(i,j)*dzt(k))  &
             *(dn0(k,i,j)+dn0(k,i,j+1))*.5
     ENDDO
  ENDDO

  k=2
  prtot=0.
  DO j=2,nyp-1
     DO i=2,nxp-1
        ppp= ( (pp(k,i,j)+pi0(k,i,j))/cp )**cpor*p00
        prtot=prtot+ppp/(dyu(i,j)*dxv(i,j))
     ENDDO
  ENDDO


  tmass=wmass+emass+smass+nmass
  aintmass=aintmass+tmass*dtlong
  area=(nxp-2)*deltax*(nyp-2)*deltay


  PRINT*,'==============================='
  PRINT*,' Mass flux - W, E, S, N'
  PRINT*,  wmass,emass,smass,nmass
  PRINT*, 'total (kg/(m2 s):',tmass/area
  PRINT*, 'total (kg/m2):',aintmass/area
  PRINT*, 'total pr change (pa):',aintmass/area*9.8
  PRINT*, 'computed mean press:',prtot/area
  PRINT*,'==============================='

  RETURN
END SUBROUTINE mass_flux
